Kuramoto model with frequency-degree correlations on complex networks 
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We study the Kuramoto model on complex networks, in which natural frequencies of phase os- 
cillators and the vertex degrees are correlated. Using the annealed network approximation and 
numerical simulations we explore a special case in which the natural frequencies of the oscillators 
and the vertex degrees are linearly coupled. We find that in uncorrelated scale-free networks with 
the degree distribution exponent 2 < 7 < 3, the model undergoes a first-order phase transition, while 
the transition becomes of the second order at 7 > 3. If 7 = 3, the phase synchronization emerges 
as a result of a hybrid phase transition that combines an abrupt emergence of synchronization, 
as in first-order phase transitions, and a critical singularity, as in second-order phase transitions. 
The critical fluctuations manifest themselves as avalanches in synchronization process. Comparing 
our analytical calculations with numerical simulations for Erdos-Renyi and scale-free networks, we 
demonstrate that the annealed network approach is accurate if the the mean degree and size of the 
network are sufficiently large. We also study analytically and numerically the Kuramoto model on 
star graphs and find that if the natural frequency of the central oscillator is sufficiently large in 
comparison to the average frequency of its neighbors, then synchronization emerges as a result of a 
first-order phase transition. This shows that oscillators sitting at hubs in a network may generate a 
discontinuous synchronization transition. 

PACS numbers: 05.45.Xt, 05.70.Fh, 64.60.aq 



T3 

C 

o 

> 
O 
OV 

in 



I. INTRODUCTION 



Synchronization phenomena attracted much attention 
of the scientific community in the last decades but the 
understanding of emergence of synchronization in com- 
plex systems is still an open problem [l[ . A few examples 
are the flashing of fireflies, the chirp of the crickets, the 
pacemaker cells of the heart, and synchronous neural ac- 
tivity. The Kuramoto model 0, y] stands out as the 
classical paradigm for studying spontaneous emergence 
of collective synchronization in complex systems (see, for 
example, Refs. This basic model is analytically 

treatable and may contribute to general understanding 
of synchronization phenomena. 

The Kuramoto model describes a system of interact- 
ing phase oscillators. An explicit solution of this model 
was found for an infinite complete graph with a symmet- 
ric single peaked distribution of natural frequencies and 
an uniform coupling constant J 0, H, 0- In this case, 
when the coupling between oscillators becomes greater 
than a critical value J c , the spontaneous synchronization 
emerges as a result of a second-order phase transition 
with the standard mean-field critical exponent j3 = 1/2 
for the order parameter. Further investigations demon- 
strated, however, that the kind of the phase transition 
depends on the form of the distribution of the natural 
frequencies of the oscillators. The Kuramoto model with 
a convex distribution function undergoes a discontinuous 
transition in contrast to the second-order transition with 
(3 = 1/2 when the distribution function is concave 
In the particular case of a flat distribution of natural fre- 
quencies, the synchronization emerges discontinuously as 



a result of the hybrid phase transition with a jump of the 
order parameter as in a first-order phase transition, but 
also with strong critical fluctuations as in a continuous 
phase transition [§-[Toj. 

Many real-world complex svstems have a structure of 
random complex networks [ill, EH, IHj], and this kind 
of structure can strongly influence their dynamics [6[. 
Within the Kuramoto model, the structure of the under- 
lying network also plays an important role and affects the 
synchronization of oscillators. The Kuramoto model with 
a symmetric single peaked distribution function of nat- 
ural frequencies on uncorrelated random scale-free com- 
plex networks with a degree distribution p(q) cx q^ 1 was 
studied in works by use of a mean-field approach. 

It was shown that if the second moment of the degree dis- 
tribution is finite in the infinite size limit (i.e., at 7 > 3), 
then the critical coupling J c is finite and the phase tran- 
sition is of the second order. In contrast to this kind of 
networks, in networks with a diverging second moment 
(i.e., at 2 < 7 < 3), the critical coupling J c tends to zero 
in the infinite size limit. This means that an arbitrary 
finite coupling leads to synchronization of phase oscilla- 
tors. A similar critical properties were found for the Ising 
and Potts models on scale- free networks 0, [13, EH • 

Recently, an interesting variation of the Kuramoto 
model was proposed by Gomez-Gardehes et al. [19f. The 
authors introduced a model in which natural frequencies 
ujj and degrees qj of vertices are rigorously (namely, lin- 
early) related, ujj — aqj + b. By use of numerical simula- 
tions of the model with N = 1000 oscillators, they found 
that a second-order phase transition occurs in Erdos- 
Renyi networks and in the model (20j . In the configu- 
rational model of scale-free networks with 7 < 3.3 and 
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in the Barabasi- Albert model they observed a first-order 
phase transition at a finite critical coupling J c in con- 
trast to the zero critical coupling found in Ref. p~34l6| 
for 7 < 3 in the infinite size limit. Gomez- Gardenes et 
al. suggested fl9| that this discontinuous transition may 
be driven by hubs that entrain and synchronize neighbor- 
ing oscillators. Recently, Leyva et al. [2l[ showed that 
the first-order phase transition is not specific for the Ku- 
ramoto model, but also can be found in other systems of 
non-linear oscillators, for example, in scale-free networks 
of interacting piecewise Rossler units. The transition was 
confirmed experimentally in electronic circuits with a star 
graph configuration 2l|. The fact that this first-order 
phase transitions can be observed experimentally opens 
ground for technological application and makes the un- 
derstanding of this behavior even more urgent. 

In the present paper, in order to understand the role 
of frequency-degree correlations for synchronization of 
phase oscillators, we carry out a detailed analysis of 
the model proposed by Gomez-Gardenes et al. [l9| in 
the case of networks with scale-free topology. Using 
the annealed network approach [(| [22j and performing 
numerical simulations of the model, we show that the 
model actually undergoes a first-order phase synchro- 
nization transition if the underlying networks have scale- 
free network topology with the degree distribution expo- 
nent 2 < 7 < 3. For scales- free networks with 7 > 3 
the system demonstrates a second-order phase transition 
with the mean-field critical exponent f3 = 1/2 for the 
order parameter. Surprisingly, we find a hybrid phase 
transition with (3 = 2/3 at 7 = 3. In the latter case, 
synchronization emerges discontinuously with increasing 
coupling between oscillators but hysteresis is absent and 
there are critical fluctuations as at second-order phase 
transitions. Interestingly, these critical phenomena are 
related to the avalanches of synchronization between os- 
cillators. Furthermore, in order to understand a role of 
hubs for synchronization, we also study the Kuramoto 
model on star graphs and find a criterion for the first- 
order synchronization transition. 



II. GENERAL EQUATIONS 

The dynamics of phase oscillators in the Kuramoto 
model is described by the following equations: 
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where N is the total number of oscillators, 6j and ujj 
are, respectively, the phase and the natural frequency of 
oscillator j, where j = 1, N. Jji > is the coupling 
between oscillators j and I. 8j is defined as 6j = d0j/dt. 
a,ji is the entry of the adjacency matrix of the network. 
ciji is equal to 1 if vertices j and / are connected, and 
dji =0 if they are not. For simplicity, we assume that 
the coupling constant is uniform, i.e., Jji = J. 



Let us use the annealed network approximation to 
solve this model on an uncorrelated random complex net- 
work [|| [22j]. Within this approach, the entries dji in 



Eq. ([TJ are replaced by the probabilities 
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that vertices j and I with degrees q.j and qi , respectively, 
are connected, (q) is the mean degree, (q) = YljQj/N- 
Here the annealed network approximation plays the role 
of a mean-field approach. Substitution of Eq. ^ to 
Eq. ([TJ means that the actual interactions of phase os- 
cillators with their nearest neighbors are replaced by 
weighed interactions with all of the oscillators. As a re- 
sult, Eq. ([TJ takes a form, 



— f2 = ojj — Q, — Jrqj siri(9j — ip), 



where 
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(3) 



(4) 



The parameter r is the order parameter of synchroniza- 
tion, tp represents a global phase of the system, and f2 is 
the group angular velocity, f2 = \&. We assume that in 
the limit t — > +00 and N — > +00, the system approaches 
a steady state with a constant group angular velocity f2, 
i.e., O = 0. 

Analyzing Eq. (j3|), one finds that there are two groups 
of phase oscillators. If — Sl| < Jrqj, then oscillator j 
is locked. In this case, Eq. ([3]) has a stable solution with 
9j = S7 and takes the form 



Q = Jrqj sin (9j — ip) 



(5) 



The locked oscillators are synchronized and are rotat- 
ing together with the same group angular velocity Q. If 
\ujj — Q\ > Jrqj , oscillator j is drifting and never reaches 
a steady state, in contrast to the locked oscillators. 

Let us study the Kuramoto model with a linear re- 
lation between natural frequencies ojj and degrees qj 
(frequency-degree correlations [HI), i.e., 

uij = aqj + b. (6) 

Using a rotating frame, u>j — > u>j — b, and rescaling the 
coupling constant, J —¥ J/\a\, one obtains that the model 
with an arbitrary parameters a and b is equivalent to the 
model with b — and a = 1. It is the case that we will 
study below. Taking into account locked and drifting 
oscillators, we write Eq. (jj) as follows, 
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where 9(x) is the Heaviside step function. The first term 
is the contribution of locked oscillators to the order pa- 
rameter and the second term is the contribution of drift- 
ing oscillators. Replacing the summation over degrees by 
integration and using an explicit solution of Eq. ([5]), we 
obtain that the contribution of the locked oscillators to 
the order parameter in the thermodynamic limit is 
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In the thermodynamic limit N — > oo, the contribution of 
the drifting oscillators to the order parameter is 
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(see Appendix lAl and @, OH)- In order to simplify our 
calculations, it is convenient to introduce a variable 



a — rJ. 



(10) 



Then, substituting Eqs. ^ and © into Eq. ([7]) and con- 
sidering the real and imaginary parts of the order param- 
eter r, we obtain a set of two equations, 

/ + oo 
dqp{q){q-ty X 

qa 



aq 
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R(a) = - r 



(11) 



(12) 



for two unknown parameters Q and a. It is convenient 
to consider fi as a function of a, — f2(a). The function 
R(a) in Eq. (jlip is defined as follows, 
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(13) 



Solving Eqs. (|TT|) and (TT2l . we find a and the group an- 
gular velocity f2. Then, from Eq. (I10[) . we find the order 
parameter r. 



III. KURAMOTO MODEL ON THE 
ERDOS-RENYI NETWORKS 




FIG. 1. Kuramoto model with frequency-degree correlations 
on Erdos-Renyi networks, (a) The function R(a) versus a 
from Eq. (|13[) for the Erdos-Renyi network with (q) = 10. 
The dashed lines display the line a/ J at different couplings 
J. The intersection of R(a) and a dashed line gives a solution 
of Eq. (|12p and, in turn, from Eq. (I10p . determines the order 
parameter r. (b) The order parameter r versus J for the 
network. 



Let us consider the Kuramoto model with frequency- 
degree correlations, u>j — qj, on the Erdos-Renyi graph 
with a given mean degree (q). In this case, the degree 
distribution is Poissonian, p(q) = {q) q e~^ / q\. The func- 
tion R{a) given by Eq. (fTB")) is represented in Fig. DJa). 
Solving numerically Eqs. (TTTT) and (TT2|) . we find that a 
non-trivial solution appears if J is greater than a critical 
coupling J c . The order parameter r as a function of ,/ 
is shown in Fig. HR). Expanding the function R(a) at 
aCl, wc find the critical behavior of r near J c , 



r cx ( J - J c )* 



(14) 



where the critical exponent /3 equals 1/2 and the critical 
coupling J c is 



J, = 



2(g) 



7Tp(ft)f2 2 



(15) 



At the critical point J = J c , the group angular velocity 
is f2 = Cl(a = 0) and can be found from the equation 



dqp(q)- 



0. 



q-Q 

At (q) 3> 1, equations (fl"5"T) and (p~6|) give 



J c = 2v / 2/ v /^y. 



(16) 



(17) 



Note that this asymptotic result has a square-root de- 
pendence on the mean degree (q) in contrast to the result 
J c = 2/[Trg(0)(q)] obtained in Ref. [l6[ for the standard 
Kuramoto model on the Erdos-Renyi networks with a 
one-peaked distribution function g(ui) of natural frequen- 
cies. 
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IV. KURAMOTO MODEL ON SCALE-FREE 
NETWORKS 

Let us consider the Kuramoto model with frequency- 
degree correlations Eq. on scale-free networks with a 
degree distribution p(q) — Ag -7 , where A is a normal- 
ization constant and Qq is the minimum degree. For this 
purpose we solve Eqs. (fTT|) and (fT2j) . The function R(a) 
given by Eq. (|13|) and our results of a numerical solution 
of Eqs. (fTTj) and (| L2[) are represented in Fig. [2] at different 
values of the degree distribution exponent 7. Note that 
in this case the function R(a) does not depend on the 
minimum degree qq. Fig. [3] displays the function Q(a) 
found from a numerical solution of Eq. (1121) at different 
values of degree exponent 7. 

Figure &b) shows that if 7 > 3 the system undergoes 
a second-order phase transition at J — J c . At J < J c , 
Eqs. (fTTj) and (fT2j) have only a trivial solution r — that 
corresponds to the intersection of R(a) and the line a/ J 
at the point a = 0. At J > J c a non-trivial solution r ^ 
emerges. The trivial and non-trivial solutions correspond 
to two intersections in Fig.^a). The solution with r ^ 
is stable while the trivial one is unstable. 

At 2 < 7 < 3 we find that the system undergoes a 
first-order transition solution at J = J c \. In the range 
J c i < J < J C 2, hysteresis takes place. In this range, there 
are three solutions of Eqs. (fTTj) and (|12|) [three intersec- 
tions between R(a) and a/ J in Fig^a)]. The non-trivial 
solution with the smallest a is always unstable. The triv- 




FIG. 2. Kuramoto model with frequency-degree correlations 
on scale-free networks, (a) R(a) versus a for scale-free net- 
works with the degree distribution exponents 7 = 3.2, 7 = 3, 
and 7 = 2.8 (dash-dot, dash-dot-dot, and solid lines from left 
to rigth, respectively). The thin solid lines represent the linear 
function a/ J for J = 1.5 and J = 2, respectively. The inter- 
section of R(a) and a/ J determines the solution of Eq. (|12[) 
and, in turn, Eq. (|10|) gives the order parameter r. Inset rep- 
resents the function <I?(a), Eq. ()19|) . for scale-free networks 
with the degree distribution exponents 7 = 3.2, 7 = 3, and 
7 = 2.8 (dash-dot, dash-dot-dot, and solid lines, from up to 
down, respectively). The thin solid horizontal lines display 
the value of 1/J for J = 1.5 and J = 2, respectively. The in- 
tersection between the solid and dashed lines determines the 
order parameter r, Eq. (|lUp . (b) The order parameter r versus 
a for scale-free networks with the degree distribution expo- 
nents 7 = 3.2, 7 = 3, and 7 = 2.8 (dash-dot, dash-dot-dot, 
and solid line from left to right, respectively). 




y=3 




y=3.2 



a 2 

FIG. 3. Group angular velocity il(a) versus a in the Ku- 
ramoto model with frequency-degree correlations on scale- free 
networks with the degree exponent 7 = 2.8, 3, and 3.2. 

ial solution a = and the solution with the largest a 
correspond to stable and metastable states (this will be 
discussed below). 

If 7 = 3, we find that there is a discontinuity in the 
order parameter r at the critical coupling but there is no 
hysteresis at J > J c . It is actually a hybrid phase tran- 
sition similar to the transition found for the fc-core of 
random graphs [13, [H[ , bootstrap percolation [29[ , and 
the avalanche collapse of interdependent networks (30j . 
Assuming that avalanches are a generic feature of hybrid 
phase transitions, we suggest that avalanche collapse of 
synchronization also occurs in the Kuramoto model at 
the critical coupling J c . When J decreases and tends to 
J c , avalanches of desynchronization of oscillators emerge 
reducing the synchronization. Namely, when a single 
oscillator becomes drifting, it triggers an avalanche in 
which a large group of previously locked oscillators be- 
come drifting. The averaged size of these avalanches 
approaches infinity as J — > J c . The structure and the 
statistics of avalanches were studied in detail for the fc- 
core problem [28[ and the collapse of interdependent net- 
works [30} . 

In order to study the case of 2 < 7 < 3, it is convenient 
to rewrite Eq. (|KT1) in the form 

Ha) = 1/J, (18) 

where we introduced the function $(a) = R(a)/a, 

$(a) = (7 - 2) (jj^) 7 J +1 dx(l - axy- z x 



v/T - ^6 (1 - ax) 8 (ax - q ° p a, J , (19) 

and the variable x = (q — fl(a)/(aq). Using Eq. (fT5|) . we 
find a criterion for a first-order phase transition. Namely, 
a first-order phase transition takes place if the function 
$(a) has a maximum at a 7^ 0. In order to prove this 
criterion, note that $(a) — > as a —¥ 00. Therefore, 
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for Eq. (fT8|) to have more than one solution, it is suf- 
ficient that 4>(a) be an increasing function of a near 
a = 0. The first derivative of $(a) is zero at a = 0, 
<£>'(« = 0) = d$(a)/(da)\ a=Q = 0, so the second deriva- 
tive, $"(a = 0) = <9 2 $(a)/(<9 2 a)| Q=0 , determines the be- 
havior of at small a. Then, the sufficient condition 
to have a first-order phase transition is <f>"(a = 0) > 0, 
that is 



( 7 -4)( 7 -3) n"(0) 



4( 7 



Q(0) 



>0, 



(20) 



where Q"(a = 0) = d 2 n(a)/(d 2 a) > 0\ a=Q . In order to 
find when the inequality is satisfied, we analyze behavior 
of fi(a) at a = for 7 close to 3. For | 7 — 3| <C 1, we 
obtain (see appendix |B1) 



n(0)/gb-2= -j(7-3), 
n"(0)/q o ~ 1.71 (7-3). 



(21) 
(22) 



Substituting these results into Eq. (|20|) . we find that 
the inequality is satisfied if 7 < 3. Solving numerically 
Eqs. (fTT|) and (|12|). we find the phase diagram shown 
in Fig. |U One can see that in the region /, there is 
no spontaneous synchronization. Spontaneous synchro- 
nization appears in region //. Region III is the region 
with hysteresis (there are one stable and one metastable 
states). 




2.5 3.0 3.5 4.0 4.5 5.0 

Y 

FIG. 4. 7 — J plane of the phase diagram of the Kuramoto 
model on scale-free networks with frequency-degree correla- 
tions. In region I (J < J c i) there is no spontaneous syn- 
chronization and the order parameter r = 0. Synchronization 
appears in region II (J > J C 2) in which the order parameter 
r > 0. Region /// (J c i < J < Jc.2) is the region of hysteresis 
with one metastable and one stable states. 

The critical behavior of the order parameter r near the 
critical point J c can be found using the Taylor series of 
the function $(a) in Eq. ([18]) at the point a = 0. We find 
that at 7 > 3, the phase transition is of the second order, 
and the order parameter r has the critical singularity (|14[) 
with the critical exponent /3 = 1/2. At 7 = 3 the model 



undergoes a hybrid phase transition with a jump r c 7^ 
of the order parameter and demonstrates the following 
critical behavior 



r - r c cx (J - J C Y 



(23) 



where the critical exponent j3 = 2/3 (see appendix [Cj) . 
The same critical exponent for the hybrid phase transi- 
tion in Kuramoto model with a flat distribution of natu- 
ral frequencies was found by Pazo 8] . This critical behav- 
ior is in contrast to = 1/2 found for hybrid transitions 
in other systems [27- 30] . Note that in the hybrid transi- 
tions, the distribution of avalanches over size S becomes 
power-law at the critical point, for example, P(S) oc S~ a 
with a = 3/2 for fc-core problem. We do not know yet 
if the exponent a takes the same value for the synchro- 
nization hybrid transition. 

Thus, the analytical consideration of the Kuramoto 
model with frequency-degree correlations on uncorrelated 
random scale-free networks shows that the type of the 
phase transition is changed at 7 = 3 from the second- 
order transition at 7 > 3 to the first-order transition at 
7 < 3. Below we will show that simulations of the model 
on the static model of scale-free networks [24]-[26| confirm 
this analytical result. This conclusion contrasts with re- 
sults of numerical simulations of Gomez-Gardenes et al. 
in Ref. in which the the first-order phase transition 
was observed even in the configuration model of scale- free 
networks with 7 s» 3.3. The reason of this disagreement 
may be related to the fact that Gomes et al. simulated 
the Kuramoto model on the top of scale-free networks 
with N = 1000 oscillators, while we simulated the Ku- 
ramoto model on networks of larger size, N = 10 4 . It 
is well-known that the clustering coefficient is finite in 
the configuration model of finite size and decreases with 
increasing size, approaching zero in the infinite size limit 
[3]], [32j- Thus, networks of small size have a larger clus- 
tering coefficient in comparison with networks of larger 
size. We suggest that clustering or degree-degree corre- 
lations in complex networks may influence the synchro- 
nization phase transition and they may be responsible 
for this discrepancy. Networks generated by the static 
model, which we use in our simulations, are uncorrelated 
and have zero clustering at 7 > 3 in the thermodynamic 
limit while weak disassortativc degree-degree correlations 
appear at 7 < 3 [26| . In general, structural correlations 
are significant for phase transitions in complex networks. 
The influence of degree-degree correlations on the per- 
colation transition in correlated networks was demon- 
strated in Ref. 1331] . 



V. COMPARISON BETWEEN THE ANNEALED 
NETWORK APPROACH AND SIMULATIONS 

In order to check the accuracy of the annealed net- 
work approach, we carried out simulations of the Ku- 
ramoto model with freq uency- degree correlations Eq. ([6]) 
for the static model |24| - |26| and compared the obtained 
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results with the numerical solution of Eqs. (fTTj) and (IT21) . 
We solved dynamical equations © by use of the Runge- 
Kutta 4th order method. In our simulations, we in- 
creased and decreased the coupling constant J and let 
the system to relax after every change of J. To find a cor- 
rect solution of Eqs. ((3]), we used a time step At = 0.001 
and a coupling step AJ = 0.02. The size of the network 
was N = 10000. 

Figure [5] displays results of our simulations for the 
Erdos-Renyi network with the mean degree (q) = 10 and 
(q) = 50. One can see that the theoretical calculations 
and the simulations agree well if the mean degree is large 
enough. For (q) = 50 the numerical results are in good 
agreement with the simulation, but for (q) = 10 there 
are some differences at J near J c . However, even at a 
small mean degree (q), (q) < 10, the annealed network 
approach give us a good description of the Kuramoto 
model. 




FIG. 5. The order parameter r versus the coupling J for the 
Kuramoto model on Erdos-Renyi networks. Numerical sim- 
ulations in the case of increasing and decreasing J are repre- 
sented by the symbols (►) and (■<), respectively. The solid 
lines represent the results of the annealed network approxi- 
mation [Eqs. (II 1 p and l|12p]. The mean degree (q) — 50 and 
size N = 10000. 

Figure [5] displays our results for a scale-free network 
generated by the static model [24| with the mean degree 
(q) = 50 and size N = 10000. As one can see, the simula- 
tions are in a good agreement with the annealed network 
approach, Eqs. (fill and (TT21 . With decreasing mean 
degree, some deviations between the simulation and the 
annealed network approximation appear, but the type of 
the phase transition is the same. Note that the critical 
coupling J c obtained by use of the annealed network ap- 
proach is slightly smaller than J c observed in simulations. 



VI. KURAMOTO MODEL ON A STAR GRAPH 

In order to reveal the role of hubs in the Kuramoto 
model on complex networks, we study the model on star 
graphs. A particular case of this system has been consid- 
ered in Ref. [l9j . In this paper, the Kuramoto model was 
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FIG. 6. Order parameter r of synchronization versus the cou- 
pling J for the Kuramoto model on scale-free networks (the 
static model) with frequency-degree correlations. Blue, red, 
and black symbols represent the results of our numerical sim- 
ulations for scale-free networks with the degree distribution 
exponent 7 = 2.9, 7 = 3, and 7 = 3.1, respectively. Nu- 
merical simulations in the case of increasing and decreasing 
J are represented by the symbols (►) and (-4), respectively. 
The solid lines represent the results of numerical solution of 
Eqs. ([TT]) and (HJ. Size N = 10 4 . 

solved explicitly in the case when the central oscillator 
has the natural frequency equal to the number K of near- 
est neighbors while the neighbors have the same natural 
frequency equal to 1. Gomez- Gardefies et al. found that 
this model undergoes the first order phase transition at 
a critical coupling [l9| . Another synchronization model, 
Stuart-Landau oscillators, was considered on a star graph 
in Ref. H|. 

Here, we obtain an exact solution of the Kuramoto 
model on a star graph with an arbitrary natural frequen- 
cies distribution in the limit of large number K of nearest 
neighbors. The dynamical equations for this model are 

K 

e,j =u>j + jY^ sin ( 9 J - & i)> ( 24 ) 
1=1 

6i=ui + Jsm(9j-0i), (25) 

where Eq. is for the central node j, and Eq. (|2"5)) 
is for its K neighbors with index I = 1,2,..., if. For 
convenience, here we define the order parameter as 

re* = i£ e *.. (26) 

Introducing this order parameter into Eqs. ([24| and (|25|). 
we obtain 

9j - =(u>j -fl)- JKr sin(0j - ip), (27) 

61 - 0) = {lui - 9) ) - Jsin(6i - 63). (28) 

The central node j is the leader and it must be locked in 
a synchronized state, 

ojj - Q = JKr sin (0j - i>) . (29) 
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where 57 = 6j is the group angular velocity. There are 
two different kinds of solutions of Eq. (|28|l for a steady 
state. If \uji — Cl\ < J , the oscillator I is locked. In this 
case, Eq. ([3]) has a stable solution, with 6i = £1 and 



— Q = Jsin(# 7 - — 



(30) 



If \uji — 0| > J, the oscillator I is drifting. Contributions 
of locked and drifting oscillators to the order parame- 
ter r can be obtained by using the method described in 
Appendix [A] Then, equation (f2"rJ)) takes a form, 



1 N 
K f-^ 



»i)Q 1 



1=1 
N 



bJl — 



Jr 



K ^ 



i=i 



LOl — Q 



Jr 



- 1 



(31) 



We introduce the distribution function of the natural fre- 
quencies of the neighbors as follows 



K 



(32) 



i=i 



Separating the imaginary and real parts in Eq. (|31[) . in 
the limit K 3> 1, we obtain a set of two equations for r 
and Q: 



( ^3 

V KJ 



dwg(ui+Q)Jl-(^j 



,(33) 



Q — uij 
K 



+ 00 



duj g (lu + x 



' ■ ■ \ / 1 \ J ~) wt.k -•/. 



(34) 



In fact, equation (|33p determines r as a function of the 
group angular velocity ft that must be found by solving 
Eq. QD. 

The analysis of Eqs. ([55)1 and shows that if the 
difference uij — (w;) between the natural frequency ajj of 
the central oscillator and the averaged natural frequency 
(lui) of its neighbors is smaller than a critical value u c , 
then synchronization between oscillators occurs at any 
nonzero coupling J. Here, (uii) = K~ x y^Li tends to 
the mean value ZJ in the limit K — > oo. Figure Eta) shows 
that, in this case, the order parameter r increases grad- 
ually with increasing J while the group angular velocity 
fi decreases. In contrast to this case, if the difference 
ujj — (uji) is larger than Li c , i.e., 



LOj - (uji) > LO c 



(35) 



then the Kuramoto model on the star graph undergoes 
a first order phase transition with hysteresis in a region 
J c i < J < J C 2- Behavior of r(J) and f2(J) is represented 
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FIG. 7. Order parameter r and the group angular velocity Q 
versus the coupling J in the Kuramoto model on a star graph. 
The results of simulations of the model for increasing and 
decreasing J are represented by blue and red lines with dots, 
respectively. We used a normal distribution of the natural 
frequencies with variance (a = 0.5) and zero mean value, 
U = 0. The number of oscillators is K = 10000. The central 
node frequency is uij = 1 < uj c ~ 2.12 for the panels (a) 
and (c), and Uj — 10 > lj c for the panels (b) and (d). The 
numerical solution of Eqs. (|33[) and (|34[) is represented by 
solid lines on these panels. 



in Fig. [7l Near the limiting points J c \ and J C 2 of the 
metastable states, i.e., at either < J/J c i - 1 < 1 or 
< 1 — J/J C 2 *C 1, r and demonstrate a universal 
critical behavior: 



r-r^cc^/Jcr-l) 1 / 2 , 

r-rc2K-(l- J/J c2 ) 1/2 , 
0-f> c ioc-(J/J c i-l) 1/2 , 
ft-ft c2 cx (1- J/J c2 ) 1/2 , 



(36) 
(37) 
(38) 
(39) 



where r c i , r C 2 , f2 c i , and Q C 2 are values of r and £1 in the 
limiting points J c \ and J C 2 (see Appendix |D|) . This first- 
order phase transition is similar to one we found in scale- 
free networks in Sees. HVl andlVl 

At uij = u> c , the order parameter r and the angular 
group velocity are continuous functions of J. However, 
at a critical coupling J c , the model has critical behavior, 



r-r c cx (J/J c - 1) 1/3 , 

n-n c <x-(j/j c -i) 1/3 , 



(40) 
(41) 

that is different from Eqs. ([3"rj ]) -([3"9" ]l (see Appendix [D). 

To check the analytical approach, we compare the nu- 
merical solution of Eqs. (f3"3"]) and with simulations 
of the Kuramoto model on the star graph. In simula- 
tions, we solved dynamical equations ([M]) and ([2"5j) by 
use of the Runge-Kutta 4th order method for the nor- 
mal distribution function g(ui) with the variance a = 0.5 
and zero mean value 57 = 0. In our simulations, we in- 
creased and decreased the coupling constant J step by 
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step and allowed the system to relax to a steady state 
after every step. We used a time step At = 0.01 and the 
coupling step A J = 0.005. The number of neighbors was 
K = 10000. Figure [7] displays our results of simulations 
and numerical solutions of Eqs. (|33j) and (f34|). One can 
see that despite some noise, the results of the numerical 
solution of Eqs. (j3"3")l and and the simulations are in 
a good agreement. 

Analyzing Eqs. (|33| and (l34l) . we find that the critical 
value lo c depends on the distribution function of natural 
frequencies of the oscillators around the central oscillator. 
In the case of the normal distribution of we obtain 
ui c /a ~ 4.25. The critical frequency uj c is calculated in 
Appendix |DJ Figure [8] displays the uij — J plane of the 
phase diagram of the Kuramoto model on a star graph. 




0-1 . 1 . 1 . 1 1 r- 

0,0 0,1 0,2 0,3 0,4 



FIG. 8. J — u)j plane of the phase diagram of the Kuramoto 
model on a star graph with 10000 neighboring oscillators in 
the case of a normal distribution of the natural frequencies 
with variance [a = 0.5) and zero mean value, u; = 0. The gray 
(orange online) region represents the region with hysteresis. 
The open dot marks the critical value uj c of the frequency ujj 
of the central oscillator. The critical frequency u c ~ 2.12 

The conclusion that the Kuramoto model on a star 
graph satisfying the condition Eq. (|35[) undergoes a first- 
order phase transition, gives a qualitative understand- 
ing of the role of hubs in the first-order phase transi- 
tion discussed in Sec. lIVI Indeed, in the Kuramoto model 
with frequency-degree correlations on a complex network, 
Eq. the natural frequency of a hub have a good 
chance to satisfy the condition (|35[). since a high degree 
of a vertex guarantees its high natural frequency. There- 
fore, if the fraction of hubs is sufficiently large, then they 
can induce a first-order synchronization phase transition. 



VII. CONCLUSIONS 

In the present paper we developed an analytical ap- 
proach based on the annealed network approximation to 
the Kuramoto model with linearly coupled natural fre- 
quencies and the degrees of vertices in complex networks 
(frequency-degree correlations). We demonstrated that 



the model undergoes a first order synchronization phase 
transition on uncorrelated scale-free networks with the 
degree distribution exponent 2 < 7 < 3, i.e., in the case 
of divergent second moment of degree distribution. A 
second-order synchronization transition occurs at 7 > 3, 
i.e., when the second moment is finite. At 7 = 3, the 
model undergoes a hybrid phase transition that combines 
a jump of the order parameter at the critical point as in 
first-order phase transitions and critical phenomena near 
the critical point as in second-order phase transitions. 
In the case of hybrid transition, avalanche collapse of 
synchronization occurs at the critical coupling J c . We 
compared our analytical calculations with numerical sim- 
ulations for Erdos-Renyi and scale-free networks of size 
N = 1000 - 10000. Our results demonstrated that the 
annealed network approach is accurate if the size of the 
network and the mean degree are sufficiently large. In 
order to understand a mechanism of the first-order syn- 
chronization phase transition, we also analyzed analyti- 
cally and numerically the Kuramoto model on star graphs 
and showed that the central oscillator plays the role of 
the leader in synchronization. If the difference between 
a natural frequency uij of the central oscillator and the 
averaged natural frequency Zo of its neighbors is smaller 
than a certain critical value w Ci i.e., OJj — 57 < u c , then 
synchronization occurs at any nonzero coupling J and it 
is gradually enhanced with increasing J. In contrast to 
this case, if L)j — uj < cj c , then the system undergoes a 
first-order transition into a synchronized state. In this 
case, hysteresis takes place in a certain range of the cou- 
pling J. This result evidences that hubs in a complex 
network of phase oscillators may play a role of driving 
force for a first-order phase transition. 
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Appendix A: Contribution of drifting oscillators to 
the order parameter 

Let us calculate the contribution of drifting oscillators 
to the order parameter r in Eq. (J7]). Despite a random 
movement of these oscillators, the total contribution of 
drifting oscillators to the order parameter becomes time 
independent in the limit t — > 00 and N —¥ 00. This con- 
tribution can be calculated by use of the density function 
p(8, qJr) that measures the density of oscillators with the 
angular velocity 9 and the mean degree q as it was shown 
in Ref. [To| . In this Appendix we present a new method 
that can also be useful for studying dynamics and relax- 
ation in the Kuramoto model. For simplicity, we calcu- 



9 



late this contribution for star graphs (see Sec. I VI I) , but with an analytic solution of Eq. (|25l) for drifting oscil- 
thc method can be generalized to other graphs. We begin lators with natural frequencies satisfying the inequality 



\uij — fi| > J. This solution is 



9/ — 6 j ■ = 2 arctan < 



J + tan 



(fc,+i)V(^-tt) 2 - J 2 /2 J( W/ -*) 2 -J 2 



UJl — Q 



(Al) 



I 

where kj are parameters determined by initial conditions Substituting Eq. (jAlj) into the last term in Eq. (|7|), we 
at given w; and J defined in the interval < kj < 2n. obtain that in the thermodynamic limit the contribution 

of all drifting oscillators is given by the equation 



J J dkduj G(k,uj + n)exp j 2i arctan + tan Vw 2 - J 2 j y/w 2 - J^jw^ 1 X 0(|w| - J), (A2) 



where we introduced the function 
1 K 



G{k,u + n) = — ^28(lj + Q - Ljj)5(k - kj). (A3) 



In order to simplify the calculations, we introduce a vari- 



able a = u}J 1 — ( J/lo) . Then Eq. (|A2[) takes the form 



"°° r+ °° „ , G(k,u(a) + Q) 

akaa exp < 2i arctan 



oo J — oo 



J + tan(i(fc + i) \a\) \a\ 



(A4) 



I 

In order to find this integral we replace the variable of variable given by a n = 2mr with n € Z, and y is a 
integration a to a = a n /t + 2y/t, where a n is a discrete continuous variable defined in the interval [0,7r]. Then 

Eq. (|A"4t takes the form 



2 
1 



E 



dkdy G(kMa Jt + 2y/t) + n) ^ ^ 



J + tan (|a n |4±* + y h±t) Q an \/t + 2y/t) 



\j l + (a n /t+2y/t) 



(A5) 



Since < y < n and < k < 2n, in the limit t — > +oo, 
we have k/t — > and y/t — > 0. Furthermore, the sum- 
mation over a n can be represented as integration over a. 
The resulting function is a periodic function of y. The 
integration over y remove the dependence on the initial 
conditions and we obtain a triple integral that does not 
depend on time and the initial conditions. Introducing 
a function g(uj + ft) = G(k,u> + Q)dk we obtain a 
well-known result for the contribution of drifting oscilla- 



tors to the order parameter r in Eq. ([7]), 



i I du g{u + (J/lo) 2 ] G(\u\ - J). (A6) 



Appendix B: Analysis of the group velocity function 

In the case of the Kuramoto model with frequency- 
degree correlations on scale-free networks with degree 
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distribution p(q) ex q~ J at q > qoi the group velocity 
function Q(a) in Eq. (fTTj) has no explicit expression. Nev- 
ertheless, it is possible to find its asymptotic behavior if 
the degree distribution exponent 7 is close to 3. In this 
case, Eq. ([TO takes the form 



(?)- ft(a) = ( 7 -l)<£ 



+ 00 

" ' / d 9 g- 7 (<?-fi(a)) x 
90 



1- 



q - fi(a) 



1- 



go 



(Bl) 



Using a new variable of integration, x = q/(q — Q), we 
can rewrite this equation as 



fi(a) 



7-2 



7-1 
9o 



_9o H(a) 

7-2 7-1 
7-3 



+00 



G?X ] X 



-3 



(B2) 



where -B = l/(51(a)/go — !)■ At 7 = 3, this equation takes 
the simple form: 



2 g0 -O(a) 
2<?o 



+00 



dx ] x 



a;- 3 v/l-(ax) 2 6 (1 - \ax\) 



(B3) 



This equation has a solution f2 = 2qo for any a. 

If 7 is close to 3, i.e., 7 = 3 + 6 where \S\ <C 1, then we 
look for a solution of Eq. (IB2[) in the form 



(B4) 
(B5) 



n(a) =q Q [2 + A] 
with |A| << 1. We find 

A = -5f a (a), 
where the function f a (ct) is defined as follows, 



/»(«)= 



1 x J 1/ 

i-Vi - o7e(i-|a|) 



(B6) 



Figure [9] displays the function / a (a). This function deter- 
mines the behavior of the group velocity function 0(a). 




FIG. 9. Function / a (a) from Eq. ([B6 



Appendix C: Critical exponent for the hybrid phase 
transition 



Let us study the critical behavior of the order parame- 
ter r of the Kuramoto model with frequency-degree corre- 
lations on scale- free networks with the degree distribution 
exponent 7 = 3. Using the solution f2 = 2qo obtained in 
Appendix IBl we find that the function $(a) in Eq. (fT9|) 
takes the form 



d>( 



a) = \ f + dxVl~x 2 9(1-| 
2 J-i 



ax\) ■ 



(CI) 



Solving Eq. (|19p , we find that a = 1 at the critical point 
,/ = J c . The region a > 1 corresponds to the synchro- 
nized state, and 



$(a) 



<ix\/l 



(C2) 



Near the critical point when 5a = a — 1 <C 1, we obtain 
<I>(a) in the leading order in <5a, 



$(i+fe)4-^(fof/ 2 



(C3) 



Expanding Eq. fjl8|> in a Taylor series in J— J c , and using 
<!>(«) from Eq. (|C2[) . we find the critical behavior of the 
order parameter r above the critical coupling J c of the 
hybrid phase transition, i.e., at J > J c , 



r - r c cx ( J - J c ) 



(C4) 



Here the critical coupling is J c = 4/7T, the jump of the 
order parameter is r c = 1/J C = 7r/4, and the critical 
exponent is /3 = 2/3. 



Appendix D: Analytical analysis of the Kuramoto 
model on a star graph 

In the case of the Kuramoto model on a star graph, we 
analyze Eq. (|34p to determine the group angular velocity 
H as a function of the coupling J. Note that the right- 
hand side of Eq. (l3~4l tends to zero with increasing K. 
Therefore, at K 3> 1, a solution fi(J) is small at small J 
and can be written as a Taylor series in J. In the leading 
order in J, equation (|34[) takes a form 



A' 

where we introduced 



A (SI) J 2 + 0(J 4 ) + -- 



(Dl) 



— 5 (w + 0) . (D2) 

-00 w 

One notes, that at the limiting points J cl and J C 2 of the 
metastable states, the first derivatives of the right- and 
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left-hand sides of Eq. (|D1[) with respect to ft becomes 
equal. This leads to an equation 

1 = KJ 2 cl(2) A' (fi cl(a) ) , (D3) 

where A'(f2) = <£4.(f2)/dO. Thus, the group angular ve- 
locities O c i and f2 c2 and the critical couplings J c \ and 
J c2 are determined by Eqs. (ID1|) and (|D3[) . Substituting 
Eq. (|D3|) into Eq. (ID1[) . we obtain an equation for f2 cl ( 2 ), 

a' (n cl(a) ) (n cl(2) - Uj ) = a (n cl(2) ) , (D4) 

From Eq. (|D3[) . we find 

J cl(2) cx 1/VF. (D5) 

Solving Eq. (|D1[) near the limiting points J c i and J c2 of 
the metastable states, i.e., at either < J/J c i — 1 <C 1 
or < 1 - J I J c2 < 1, we find 

fi(J) =fi c i- 5i(J/J c i - 1) 1/2 , (D6) 
fi(J) =O c2 + S 2 (l- J/J c2 ) 1/2 , (D7) 



where B 1(2) = 2|A(fi cl ( 2 ))/-A" (fl c i(2)) | 1/2 - According 
to Eq. (|3"5j) . the order parameter r also has this kind of 
singular behavior near J c \ and J c2 . 

Now let us find the critical frequency w c of the central 
oscillator at which hysteresis disappears, i.e., J c \ — J c2 = 
J c . Analyzing Eq. (|D1[) . we find that J c , O c , and tu c can 
be found from a set of equations 

n c -u c = KJ 2 C A (Q c ) , 

i = k j 2 c a' (rg , (D8) 
a" (n c ) = o. 

For the gaussian distribution with variance a and zero 
mean value we find uj c /a ~ 4.25. In order to find critical 
behavior near this special point, < |1 — J/J c \ <C 1, we 
solve Eq. (|D1[) and find 

Q{J) = n c -B 3 (J/J cl -l) 1 / 3 , (D9) 

where B 3 = \12A(Cl c )/A'" (O c ) I 1 / 3 . Thus, at Wj = w c , 
the order parameter r and the angular group velocity are 
a continuous function of J but, at J = J c , they have a 
singular behavior Eq. (|D9[) . 
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